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Abstract 

Certain types of electro-magnetic waves propagating in a plasma can undergo a mode conversion process. In magnetic 
confinement fusion, this phenomenon is very useful to heat the plasma, since it permits to transfer the heat at or near the 
plasma center. This work focuses on a mathematical model of wave propagation around the mode conversion region, from 
both theoretical and numerical points of view. It aims at developing, for a well-posed equation, specific basis functions to 
study a wave mode conversion process. These basis functions, called generalized plane waves, are intrinsically based on 
variable coefficients. As such, they are particularly adapted to the mode conversion problem. The design of generalized 
plane waves for the proposed model is described in detail. Their implementation within a discontinuous Galerkin method 
then provides numerical simulations of the process. These first 2D simulations for this model agree with qualitative aspects 
studied in previous works. 

Keywords: Wave propagation, variable coefficient, mode conversion, generalized plane waves. 


Introduction 


Q Mode conversion corresponds to a transfer of energy between different types of propagating waves. It is an important 
• ^ problem in magnetic confinement fusion particularly for plasma heating or current drive. Indeed, some waves used for these 
misapplications cannot propagate directly from the launching region at the wall toward the point of the plasma where they 
would be useful. But some other forms of wave can be sent from the wall toward the plasma, to be converted into the desired 
Q^wave at the mode conversion region, and so penetrate until the heating or control point. See [20] for a first study of mode 
conversion equations. Even though experimental models [angnaiiH] and simple one dimensional models 0122] have been 
(N studied, the two dimensional mathematical model is still not well understood. In [24] a two dimensional model is derived by 
^ means of an asymptotic expansion, but the resulting equation is not standard in the literature ; propagating solutions are 
then constructed thanks to an integral representation. 

Mode conversion corresponds to a propagating wave transmitted from one propagative zone to another one, even though 
the two zones only touch along a curve. The mode conversion region is defined as the vicinity of this curve. The two 
dimensional model studied in the present work comes from the cold plasma model for wave propagation in a plasma confined 
by a magnetic field, and reads 

o 


{dl + {d + d)d^dy + \d\^dl) F+{d- d)xdyF 


IH-h x{x -h y) 

jl 


F = 0, 


( 1 ) 


. ^ where F is the scalar unknown and d and /i are complex parameters linked to the confining magnetic field and the electron 
density. Most of the derivation will follow the steps of [24]. However a different exposition is proposed here, highlighting the 
?-H different steps of the reasoning and insisting on a crucial change of unknowns. Moreover, Equation .§ does not appear in ED, 
^ where an equivalent equation is given for a different unknown. Compared to the latter, Equation 0 has a particular structure, 
more convenient to prove the well-posedness. This work presents the first well-posedness result for a mode conversion equation 
obtained by expansion in the mode conversion region. In this elliptic second order linear 2D partial differential equation, 
the physical properties of the domain appear in the zeroth order term, as the sign of f{x^y) = x{x y). The domain is 
propagative on {(x^y)/f{x^y) < 0} and absorbing on {{x^y)/f{x^y) > 0}, see Eigurej^and the mode conversion region is 
reduced to the vicinity of the origin. It is clear that mode conversion is strongly linked to variable coefficients. 

This work aims at developing Generalized Plane Waves (GPWs) to study a wave mode conversion process. GPWs have 
been developed following the idea that plane waves are relevant basis functions to solve wave problems numerically, since the 
oscillatory behavior of the problem is embedded in their definition. In the same way, GPWs encode information about the 
problem to be solved, but they are specifically adapted to problems with variable coefficients. They were introduced in m, 
and their interpolation properties for Helmholtz equation with a variable coefficient were presented in m- In this work. 
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Figure 1: Level curves of the function /(x, y) = x{x + y)^ indicating the propagative (P) and evanescent zones as well as the 
two cut-offs lines, x = 0 and x + ^ = 0, limiting the different zones. They cross at the mode conversion point (0, 0). 


for the first time, GPWs are designed for a mode conversion equation, namely 0’ and are implemented in a discontinuous 
Galerkin method for numerical simulation of the mode conversion process in this model. 

First numerical evidences of a mode conversion process are displayed for this 2D full-wave expansion model. A typical 
test case is proposed together with a way to estimate the transmission coefficient across the mode conversion region. The 
influence of different parameters is illustrated, mainly the influence of the incident angle already studied in [20], in a ID 
model in [22] and in the 2D propagating solutions in [24] . 

Section presents the derivation of the second order equation 0, describing the mode conversion process, while the 
existence and uniqueness of a weak solution to this equation are proved in Section [^ Even though Section [^ is self contained 
and does not require any prior knowledge on waves in plasma, it is independent of the rest of the article. Therefore it can 
be skipped by a reader willing to start from Equation 0. This work then focuses on numerical aspects. Section first 
describes a discontinuous Galerkin method for numerical simulation. Then GPWs adapted to the mode conversion equation 
are carefully designed, distinguishing between the individual construction of an approximated solution to the equation and 
the global features of a set of independent GPWs. Numerical results are finally displayed in Section [^ showing evidences of 
mode conversion and highlighting the influence of different parameters on the process. 


2 A wave propagation model in the mode conversion region 


This section presents the formal process leading to the equation studied in the rest of this article. It is mainly based on 


[24] . since Equation (22) is Equation (59) from that reference. The goal of this presentation is to shed a different light on 
the derivation of the equation. Eor the sake of completeness, preliminary material related to plasma physics and the wave 
propagation model is presented first. The following subsection then focuses on Maxwell’s equations in the mode conversion 
region: since the dielectric tensor is defined in a simpler way in a specific orthogonal coordinate system, the idea is to obtain 
a reduced system with fewer components of the electric and magnetic fields in those coordinates. Such a simplification follows 
naturally from the techniques of geometrical optics, see m- It is based on physically relevant hypotheses concerning the 
relative order of the studied quantities, and subsequent expansions valid in the mode conversion region. Gonsider e, the 
inverse of the classical geometrical optics expansion parameter. If uj is the wave frequency, c the speed of light in vacuum 
and L the characteristic equilibrium wavelength, e is defined as e = cI{ujL). The model will result from a formal expansion 
with respect to e. The scaling assumptions follow the perpendicular stratified case, as studied in [25] . 


2.1 Preliminaries 

The toroidal geometry can be described by the classical axisymmetric coordinates (r, 0, z) G x [0, 27r) x M, see Eigurej^ 
and the corresponding orthonormal right handed basis {er,e(f),ez). The poloidal plane is defined as a half plane given by 
a constant 0, so that {er^ez) is an orthonormal basis of the poloidal plane while is orthogonal to the poloidal plane, 
and = rV(/) defines the toroidal direction. However, since the wave propagation phenomena are driven by the confining 
magnetic field, a toroidal coordinate system adapted to the magnetic field will be useful to describe the wave propagation 
model. We assume that the magnetic field is known from an independent solution of MHD equilibrium. The magnetic field 
lines display a helical shape, winding around the interior of the torus and that the flux surfaces are closed and nested. In 
such a case, flux coordinates form a set of coordinates adapted to the shape of the flux surfaces of the confining magnetic 
field. They consist of a radius-like coordinate, the flux label ip G ['0min, V^max] , and two angle-like coordinates, the toroidal 
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Figure 2: Left: Axisymmetric coordinates and toroidal geometry. Right: Poloidal angle 0 and level curves of the flux label 
7/^ in a poloidal half plane. In the poloidal plane, the magnetic axis is the point enclosed by all the flux surfaces. The mode 
conversion point is distinct from the magnetic axis. 


angle 0 G [0, 27r) and poloidal angle 0 G [0,27r). The bounds on -0, '0min and V^max, respectively correspond to the magnetic 
axis and the outermost closed flux surface. The toroidal angle (j) is the same as the axisymmetric coordinates angle, and 
the coordinates in the poloidal plane are described in Figure The magnetic flux labeling a curve in the poloidal plane 
measures the flux of magnetic field across the surface enclosed by the curve. It is increasing from the magnetic axis to the 
boundary of the plasma, and is such that is orthogonal to the confining magnetic field b. The poloidal angle 0 is such 
that (VV^, V6>) is a non-orthogonal basis of the poloidal plane, with the same orientation as (e^, e^). The resulting covariant 
basis (VV^, V^, V0) is a left handed non-orthogonal non-normalized basis, associated with the flux coordinates ('0,^,0). 

The mode conversion region, as the intersection of two cut-off surfaces, has been introduced as a curve. It intersects 
each poloidal plane at a single point, the point (r, = {vq^zo) on Figure]^ which stands off the magnetic axis. The mode 

conversion point is the origin of the flux coordinate system in the poloidal plane,that is to say ('0(ro, 2 : 0 ), 6>(ro, 2 : 0 )) = (0,0). 
Since the origin is not the magnetic axis, at the origin ('0,6>) = (0,0) there is no problem to define VpJ while V6> is not 
defined. The model will be derived based on a formal expansion with respect to the small parameter e, the inverse of the 
classical geometrical optics expansion parameter. Following [24] , the mode conversion region is characterized by variations of 
the poloidal coordinates ( 7 /^, 0) of order yT, while the variations of the toroidal angle (p are of order 1/e as well as er = 0(1). 
In this regime the poloidal coordinates can be written as ( 7 /, 0) = ('0(er, ez)^0{er^ ez)) where pj and 0 have derivatives of order 
0(1), so that the derivatives with respect to the flux coordinates will be scaled as e(9y,90,90), and likewise the gradient 
vectors will be scaled as ( —, ^ ^). 

The curl operator is crucial in electromagnetics. In order to express it in the flux coordinate system in a concise way, 
consider the scaled contravariant basis associated with the ('0,^,0) coordinates, defined by 

(VO Vp Vp Vip Vip V0\ 

Ve e e e e e J 


The confining magnetic field b has a poloidal bp and a toroidal bt components, such that b = bp + bt- They are defined by 


V^p V(p ^Vp; VO Vtp V(p _ VO 

b = - X - Q - X — , where bp = - x -, b^ = —Q - x —, 

eeee e e e e 

and the safety factor only depends on the magnetic flux, that is to say Q = Q{p;]^FoT clarity 6, bt and bp will respectively 
denote |b|, jb^l and |bp|. The Jacobian of the covariant basis thus reads ^ x ^ = ~^t/(^^Q)- The contravariant basis 


^The safety factor measures the winding of the magnetic field lines around the torus. It is proportional to the ratio between the toroidal and 
poloidal fields The word safety refers to the resulting stability of a configuration, since a high Q tends to improve the stability and therefore 

the safety. 
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also provides an exact expression of the curl operator: 


V X V = e X ^{deV^ - d^Ve) + ^ x - d^V^) 

+~ X ^{d^Ve - deV^)^ . 

Another basis will be used to give a simple expression of the dielectric tensor. This other basis is right handed and orthonormal, 
and is defined by a first vector aligned with while the third vector is aligned with the magnetic field b: 

, 02 = b X ei/ 6 , and es = ei x 02 = h/h. 


ei = 






The different components of any vector V will be denoted as 


VO Vd) 

V = = Viei 

e e e 


■ V2e2 + Vs^s- 


( 2 ) 


Because the covariant, contravariant and orthonormal bases depend on the space variables, any derivative of a vector 
quantity expressed in any of these bases will involve a derivative of the basis vectors, which becomes a remainder of order e 
thanks to the scaling in the mode conversion region. For example one has 

ed^ + Ved^^ = eO^V^ + 0 (e), 


so that in the covariant basis the divergence operator reads 




V 6 » 


V0. 


v-v =e{^d^ + —d0 + ^ds]-{—v^ + —V0 + ^v, 


e e e 

= S5iVi + © 2^2 + © 3^3 + 0 (e), 




V 6 », 


V(/., 


(3) 


where the differential operators 2 ) are defined by 


2)i 


ei • V = eei • 




, ©3 


63 • V 


bt 

erQb 


e{d0 + Qd^), 


2)2 = ©2 • V 



e {-b^de + Qbld^,) . 


Since the equilibrium is assumed to be axisymmetric, the fields vary as where N satisfies eN = 0(1). This essentially 
leads to a 2D reduction of the 3D model by restricting the study to the poloidal plane. In this context the mode conversion 
region then becomes a point, as the intersection of two cut-off curves, and lies at the origin of the poloidal plane (-0, 0) = (0, 0). 


2.2 The cold plasma model 

We consider here a toroidally confined plasma. In an axisymmetric equilibrium state we study an incoming wave propagating 
in the poloidal plane as a linear perturbation, reducing the model to a 2D problem. Different mode conversion processes exist, 
and this work focuses on mode conversion between the so-called ordinary (O) and extraordinary (X) propagation modes. 
These propagation modes can be described in terms of components of the wave electric field with respect to the direction of 
the confining magnetic field : a pure 0-mode wave electric field only has a parallel component, while a pure X-mode wave 
electric field only has perpendicular components. The X-mode wave considered in this work is a left-handed polarized wave. 

The cold plasma model corresponds to propagation of time harmonic electromagnetic waves through zero-temperature 
plasma. Maxwell’s equations are combined with a linearized momentum equation for the particle motion in a stationary 
confining magnetic field b. The thermal speed is neglected with respect to the wave speed. Even though this work is focused 
on high frequency waves, this model encompasses a much broader range of wave motion than magneto-hydrodynamic models. 
The coupling between the electro-magnetic fields and the fluid motion appears via the current generated by particle motion, 
modeled as a source term in Maxwell’s equations. Frequencies will be expressed in uj units while distances will be expressed 
in c/uj units, where c stands for the speed of light in vacuum. The resulting time-harmonic system reads 


V X E = -^B, (4) 

V X B = tu'E, (5) 
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where k is the dielectric tensor for the cold plasma model. It is classically expressed in a right handed orthonormal basis 
whose third vector is aligned with the magnetic field b, so in particular in the (ei, 62, 63) basis: 


K = 



— IK.A 0 \ 

0 I 

0 K||,/ 


(6) 


where all the coefficients are varying in space. 

Following the analysis of the dispersion relation for the cold plasma model, we can describe different cut-offs as surfaces 
between propagative and evanescent zones for a given type of wave. The corresponding type of wave, impinging from the 
propagative zone on a cut-off, would then be reflected by the cut-off toward the propagative zone. The 0-mode cut-off is the 
surface defined by the condition =0 while the left X-mode cut-off is the surface defined by the condition — Ka =0. 
A pure 0-mode wave can only propagate if >0 while a pure left-handed polarized X-mode wave can only propagate if 
> 0. The mode conversion occurs in a neighborhood of the intersecting curve of these two surfaces. The model 
derived in this work was developed in [24] , and relies on an expansion in the mode conversion region. The resulting equation, 
namely inherits from the cold plasma model the fact that it has variable coefficients. 


2.3 A differential system in the mode conversion region 

This paragraph describes the reduction of the 6 x 6 system describing Maxwell’s equations 0 and 0 to a 2 X 2 system, by 
eliminating some convenient unknowns. Moreover a further simplification is performed by specifying the phase of the desired 
solutions. 

The first idea is to eliminate the (F^i, £^ 2 ) components. To that purpose, combining the 62 and es components of Faraday’s 
law one can show that these two components satisfy 

fVU^ fVU 

± ^E2) = (2>1 ± t^2)E^ T {Bi 

As a result, the first two components of 0 yield a pair of equations independent of the Ei and E 2 components: 


T ^a) 


(Di ± T (5i A ^^ 2 ) + ^Bs 

= ± iB2) + ~ '^^2)^3 + 0 (e). 


( 7 ) 


To complement the latter into the system that will later be expanded in the mode conversion region, express the third 
components of Faraday’s and Ampere’s laws in the orthonormal basis, as well as the divergence free condition for the 
magnetic field: 

d'lpE^ = —iBi ^^ -\-T)sE(p^ ( 8 ) 

ihz\\Es =—'D 2 B 1 E'D 1 B 2 0{e), (9) 

+ 1 ) 2^2 + + 0 (e) = 0 . ( 10 ) 


Equations 0 ,0,0 and ( [1Q| ) form a simplified system in which the toroidal component of the electric field, still 
appears. A series of hypotheses will then allow further simplification of this system, by scanning the relative orders of the 
different terms to identify the leading order terms. 


Hypothesis 2.1. The wave amplitudes of the fields E and B are expeeted to vary faster than the equilibrium seale length, 
hut slower than the wave variation f^/e, 0/e. More speeifieally : 


even if this representation of the magnetie and eleetrie fields is not unique. 

Definition 2.1. The phase will be denoted : ^ loeal wave number is then defined as 


k = eVA = eV 




Taking the derivative of any eomponent o/k with respeet to (j) then reduees to a simple multiplieation by eN. 
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Hypothesis 2.2. The components of the local wave number that are perpendicular to b in the mode conversion region satisfy 


ki = 0(\/e) and k 2 = 0 (\/e), 

( 11 ) 

while the parallel component ks = k\\ satisfies 


;„2 _ 1 ‘^pe 

" 1-fice’ 

( 12 ) 

which corresponds to the usual X-mode cut-off condition 


K± - K/\ - kl = O(v^). 

(13) 


Moreover the phase X is of order 1. 


Hypothesis 2.3. The components E(fy, Bi and B 2 scale as 0(1). 

Thanks to the order of the two first components of the wave number, see Equation 0 from Hypothesis |2.2[ then 
ks = + 0{^/e). As a result and D 2 are operators of order 0{^/e), while ^Ds is of order 0(1). Then since Equations 

(10) and 0 show that {n_\_ ^ — kl){Bi ± 2 ^ 2 ) = O(yT), Equation ( [T^ from Hypothesis |2.2| provides a way to eliminate 

the B 2 component since: 

B 2 = —iBi + 0(\/e). 

Two more unknowns can be eliminated thanks to Equations 0 and ( P!q| ): 

erh 

= eNEs + + O(V^), 

tk^B^ = —T)2B2 — + 0 (e). 

At this point, H 2 , and Bs are expressed explicitly in terms of the two last unknowns, namely Es and Hi, while Es 
and Hi satisfy an 2 x 2 differential system independent of the other unknowns. A new unknown is finally defined in order to 
obtain a diagonal differential operator : 

Definition 2.2. Define the new variable H 3 = H 3 + 2erbp/{beN)Bi. 

As a result, we now get the following system: 

erbry 


(J)i — i^2)Ei — —/^||H 3 + 2/^11 ^^^Hi + 0 (e), 

~ erb ~ 

(Di + i'D2)Es = —2/^11 ^^^3 


(14) 


2 Hi + 0 (e). 


The key point of the next step of the simplification process is to specify the phase function of the electric and magnetic 
fields introduced in Definition |2.1[ in order to cancel some of the differential operators terms and obtain a simpler system on 
the amplitude of the unknowns H 3 and Hi. Notice that specifying the phase function actually means looking for only a class 


of solutions to System (14) 


In order to simplify the T )2 term, for any component A of the electric or magnetic field. A' standing for the amplitude of 
A. Erom the definition of the differential operator, one has 


D2A = 




QK 


Qb 




- (v^dffA' + zdffXA') + ^zeNA'J 


SO that, choosing dgX = 


eNQbl 


Definition 


2.1 


, it becomes D 2 A = 

(0,0) 


VedeA 




0{^/e). Indeed the scaling introduced in 


^ensures that ^/cdeA' = 0(1) while the remaining terms scale as = o(l). Eollowing the same idea, 

the equivalent simplification of Di is obtained by setting ei • ^d^pX + ei • ^OqX = 0. Choosing the corresponding value 

of dypX, it yields DiA = y^ei • (^^d^pA' + ^deA'^ e^ + o(l). To summarize: 

Definition 2.3. From now on, the solutions of (p3 are more specifically sought with 


X{f;,0)=X{0,0) 


eNQbl 




a I 

e - ei- 


( 0 . 0 ) 




-^eNQ^\ 


b^ 




(0,0) 


6 
























At last the definitions of the O and X mode cut-offs in the mode conversion model lead to the final simplification idea: 


Definition 2.4. The 0 and X-mode cut-offs are respectively defined by 

= 0{^/e) and [n± - /yfe.O/^/e) = 0{y/e), 

see Hypothesis \2.^ The mode conversion point is uniquely defined as the point where = 0. 

Notice that each of the coefficients appearing in the right hand side of System (14) is proportional to a quantity vanishing 
at one cut-off: either or n_\_ — — k^. A spontaneous change of variables and the subsequent differential operators are 

then: 


Definition 2.5. The constant ci is defined by the local expansion - 

2 (Sn) ^11 + - «A - A:|) = Ve {C 2 tp/V^ + C 36 '/Ve)- In other words 


y/ecifi/y/e, and C 2 and cs are defined by 


Cl 


-2d^n\\{0) 

Te ’ 


a^F(o,o) a0F(o,o) 

C2 — - T - 5 C3 - 






where F{ 


jL _t) =2 


v/i’ Vi 


\bcN J 




Define the change of variables 


y'= {0 + C 2 ip I C3)e 
x' = 


while a = ^ and the constants are d'^, o-nd d'^ defined by the related operators 


= ei • f — + —i dy' ] = d[dx^ + 



V e C3 


e / 


dy, =dsdy>. 


( 0 , 0 ) 


Thanks to these definitions, the leading order terms of System (14) give the following system 

f (D'l + ^Tl2)E'^ = cix'aE'^ + C3y'{2B[), 

1 (D;-^2)^)(2Bi) =cix'E!^-cix'a{2B[), 

where it is crucial that the unknowns are now the amplitudes of the original unknowns. For the sake of clarity, another 
rescaling change of variables is introduced by the following definition. 

Definition 2 . 6 . Define ei and 62 such that ei = a/-^ and 62 = Define a scaling of variables and unknowns by 


aci 


X' = cix and y' = 62 ^, 

\ 2aB'i = B and E'^ = E, 

together with d 2 = ; d^ = and d = d 2 tds . 

The resulting system finally reads: 

{dxEddy)E =xETyB, 

{dxddy)B =xE — xB. 

It is a 2 X 2 differential system model in the mode conversion region. An important feature of this system is that, in the (x, y) 
coordinates, the cut-offs are defined by x = 0 and x y = 0. It is a direct consequence of the fact that in the intermediate 
variables they respectively correspond to x' = 0 and cia‘^x' + c^y' = 0, combined with the final rescaling. As to the mode 
conversion region, it is the neighborhood of the point {x = 0}n{x-h^ = 0}, that is to say the origin {x^y) = (0, 0). 

Since from now on the only remaining parameter is d G C, define df = 2s{d) and dr = 5R(d). 


(15) 


2.4 A second order equation 


In order to write System (15) as a single second order equation, one would naturally eliminate one of the components and 
obtain one of the following equations: 

{dx + ddy + x) (-{dx + d5j^)E ) = {dx + ddy + x) (-E ] + xE, 
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{dx + ddy - x) y-{dx + ddy)Bj = - {d^ + ddy - x)B. 

However both of these equations have an artificial singularity at the mode conversion point. Obtaining a well-behaved second 
order equation for System then requires further investigation. 

Looking for a physical solution, one can perform a change of variables involving the exponential of a quadratic form, 
as proposed in [24] (Equation (54)). Suppose that a quadratic form Q{x,y) = ^{Kx‘^ -1- 2Lxy + My‘^) satisfies (E, B) = 
(E, The idea is to determine Q such as to simplify the differential system satisfied by the amplitude functions (E, B) 

and such that 5RQ(x,^) < 0 to guarantee a physical behavior at infinity. 

Starting from System (15), it is straightforward to see that the amplitudes (E, B) satisfy 


r {d,^ddy)E ={{l-K-dL)x-{L^dM)y)E^yB, 

\ (5^ + ddy)B = xE - ((1 + + dL)x + dM)y)B. 

A simplification of the right hand side is then performed thanks to an adequate set of constants (LC, L, M), setting 

L T dM = 1 -f- K T dL = 0, 
l/y = K^dL-l = 1/{L + dM), 


so that System (16) becomes 


{dj, + ddy){yE) = -xE + yyB, 
(dxddy)B =xE — yyB. 


From equations 0 and ( [T^ , it is straightforward that d/r = 2 + l//i. 


K= 


ddy 


d — d 


L = -^ 


dy 


and M = — 




So the quadratic form reads Q{x^y) = | 


1 I ddfi \ 

d-d) ' 


d — d 




d — d 


(16) 

(17) 

(18) 

(19) 

( 20 ) 


infinity is linked to the sign of J?Q, since if 3?Q < 0 there is no propagation for |(x, y)\ 
to 


The physical behavior of the solution at 
oo. A little more algebra then leads 


23?Q(.,.) = -|| 




1 




{\iJ,fy + xy 
y 


The parameter y satisfies the second degree polynomial equation dy^ — 2y — 1 = 0. The two roots of the polynomial 
X‘^ + 2X — d are {1/y)± = —1± Vl + d, and the inequality 5RQ < 0 is satisfied if and only if di and have opposite signs. 

Since the sign of ^Vl + d is the sign of d^, it is then clear that the root (l//i)_ = —1 — Vl + d ensures the desired estimate 
at infinity for 5RQ. The following lemma summarizes this result. 


Lemma 2.1. Given the eomplex numbers {K±^ L±^ M±) satisfying (20) with 1/y eomputed as {l/y)±, and the eorresponding 
quadratie forms Q±{x^y) = ^{K±x‘^ + 2L±xy + M±y‘^). The form Q_ eorresponding to the root (l//i)_ = —1 — ^/r+~d 
satisfies 

^Q-{x,y)<0, V(x,y)GM^. (21) 

As a result the funetion expQ_(x,^) is bounded. 


So define now (E, B) = (E, B) exp(—Q_(x, ^)). Since the two right hand sides of System (19) are the same up to a 
multiplicative constant, a divergence free condition holds : dx{yE + B) + dy{ydE + dB) = 0. As a result there is a potential 
ip such that —dyip = yE B and dxT = ydE + dB. For the sake of simplicity, define the differential operator ;D = 9^, + ddy. 


Then F = 




fi(d—d) 


satisfies E = DF, B = —y'DF and the second order differential equation 


DDF+ -DF- 
F 


yyT2)F = 0. 


( 22 ) 


At this point this work finally diverges from [24|. Since we wanted a well-behaved second order equation for System (iT^, 
we can undo the transformation, to get an equation for F = FexpQ_(x,^). From ( [^ stems that DQ = ^1 -h ^ 


X, 


DQ = yy—x and DDQ = dy—1 = 1 + ^- The definition of F implies that DF = (DF—FDQ)e ^ and DF = (DF—FDQ)e 
Then (22) directly gives 


DDF -f 2idixdyF — 


IH - h x{x + ^) ) F = 0, 

y 


which is precisely Equation 0. 
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3 Theoretical study 


In order to work on a well-posed problem, a weak formulation will be derived focusing on coercivity and symmetry properties. 
Remember that d = + idi is a complex constant and that l//i = — 1 — Vl + d. 

Consider the second order term of the elliptic equation 0- The corresponding differential operator can be written in a 
divergence form, as 


2)D = {d^ + ddy){d^ + ddy) =dl + {d+ d)d^dy + \d\'^ dl = V 


1 d 
d |(i|^ 


V 


(23) 


Since the eigenvalues of the matrix M = ( i |^2 ) 0 1 + |d|^, the lower bound estimate associated with the form 


(23) reads: VX G C^, MX • X > min(0,1 + |dp)||X|p. But because min(0,1 + |dp) = 0 this is not an ellipticity condition 
for the D'D operator. However this operator can also be written 


DD = {d^ + ddy){d^ + ddy) 2drd^dy + \d\^ = V ' 

1 


1 dr 

dr |dp 


V 


(24) 


dr 

dr |dp 


are 


Since the eigenvalues of the matrix A = 

^ l + |d|2±y(l-|d|2)2 + 4d2 l + |d|2±y'(l + |d|2)2_4rf2 

A+ = --- = ---. 


the lower bound estimate associated with the form (24) reads: VX E C^, AX • X > A_||X|p. Consequently it is sufficient to 
suppose that d^ 7 ^ 0 for the latter to be an ellipticity condition for the D;D operator. So the weak formulation will be based 
on the identity 'D'Du = V • {AXu). 

As to the first order term of 0 it will be split in the weak formulation in a symmetric way, thanks to the fact that 
it can be written as 2idi\' • VF with v = I ). Indeed, this vector field v is such that the x component does not depend 


on the X variable and the y component does not depend on the y variable. So this first order term can be symmetrized as 
/^(v • VF)Tp = I /q(v • VF)Tp -h ^ Jq F(v • V^). Moreover, this justifies the choice of Equation 0 over ( [ 2 ^ since in the 
latter the first order terms ^ and yy'DF • ip do not have such a simple symmetric formulation. According to 

these ideas, complement Equation 0 with an appropriate boundary condition: 


V • (AXu) + 2idixdyU — ^ + x{x + y^ = 0, (f2), 

v • (AXu) + idiXVyU + icru = {dO)^ 


(25) 


being a Lipschitz domain around the mode conversion point, i.e. (0,0) G O, a being positive constant, and 1/y = 
-1 - A/TTd. 

Definition 3.1. Consider a Lipschitz domain Q around the mode conversion point, i.e. (0,0) E Ll. A weak solution of the 
partial differential system (25) is: for all v E 


/ AXu • Xv — idi j X {dyuv — udyv) + / c(x, y)uv Ficr uv= gv, 

Jn Jn Jn Jon Jon 


(26) 


where u E F[^{LL) and c{x, y) = l-\- pF x{x + y). 


Under a single hypothesis on the imaginary part of the parameter d, the following theorem states the well-posedness of 


(1^). 

Theorem 3.1. Consider a Lipschitz domain ft. Suppose di < 0 and a > 0 ; then there exists a unique solution u E H^{ft) 
to the weak formulation (26). 

Proof. The idea is to first deal with the first order term and then treat the whole problem as a coercive plus compact 
decomposition. 

Define the intermediate problem 


/ 

Jn 


AXu • Xv ■ 


idi X {dyUv — udyv) F \ I uv F iG I uv 

Jn Jn Jan 

= 9 V+ fv, 

Jan Jn 


( 27 ) 
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where the parameter A is to be tuned to ensure the well-posedness of this problem. Define 


It is then classical to write 


i(u^ v)= AVu • Vn — idi x {dyuv — udyv) + A / uv^ia uv. 
Jn Jn Jn Jan 

u)= AVu • Vu -i-di x'A (dyuu) + A / 

Jq Jq Jq 

>X-\\Wuf-\di\ 

and setting A = > (M-l 


> A_||Vnf - Md JxIlIVt^lllMI + A|H|^ 
{x,y)en 


( 28 ) 


(29) 


• considering A_ - |di| niax(^,y)£n k| + A, 

then 'Sia{u^u) > f 

• considering A (|p^) - |di| niax(^,y)gn + A_, 

then ?R:a{u,u) > ^||Vi4|p. 


As a result ^a{u^u) > there exists a unique solution to (27) thanks to Lax-Milgram theorem. 

Define T : (f^g) G L‘^{Q) x L‘^{dQ) u ^ u being the unique solution to ( [27| ). The operator T is compact since the 


unique solution rt of (27) is actually in and the embedding ^ L^{Q) is compact. 

Then u is solution to the initial problem ( [^ if and only if u satisfies u = T((A + c{x^y))u^g)^ which is equivalent to u 
satisfying (/ — T((A + c(x, ^))-, 0))u = T(0, g). Note that the composition function T(', 0) composed with the multiplication 
by the smooth function c(x, g) + A is compact, as the composition of a smooth and a continuous functions. Since this is a 
coercive plus compact decomposition, the Fredholm alternative states the equivalence between existence and uniqueness of 
a solution. To prove the uniqueness start from a solution u to the homogeneous equation (26), i.e. with ^ = 0. Then the 


imaginary part of (26) reads 


4 / +'^ [ 

Jn JdQ 


Since a > 0 , 


< 0 and di 

fj, 


kl" =0. 


< 0 , it implies that u = 0 . So the solution is unique. 


(30) 


□ 


Note that 7 ^ 0 is crucial to prove the coercivity of the bilinear form, while the sign of di is crucial to prove the 
uniqueness. 


4 A numerical method 

The main feature of System ( |2^ is that two of the coefficients in the equation depend smoothly on the space variables. 
Moreover the mode conversion point itself has been defined as the point satisfying both the X- and 0-mode cut-off conditions, 
namely x y = 0 and x = 0. Since each cut-off is defined as the level curve of a smooth function, it is then clear that the 
varying nature of these quantities is crucial to the mode conversion phenomenon. As a result it is important to use a numerical 
method adapted to variable coefficients. 

In order to reduce the pollution effect, documented in [T], appearing in finite elements methods for wave propagation, 
plane wave methods use basis functions adapted to this particular application: these basis functions are solutions of the 
homogeneous equation, [9] . See m for the first description of these methods under the denomination Trefftz-based methods, 
and lais] for more recent developments. The leading idea is that basis functions embedding information about the problem 
of interest are more efficient than polynomial functions to approximate a wave. 

The numerical method that we propose in this work relies on basis functions designed to fit the variable coefficients, called 
Generalized Plane Waves (GPWs), coupled with a Discontinuous Galerkin (DG) method, called the Ultra-Weak Variational 
Formulation (UWVF). The UWVF, proposed by B. Despres in [8], has been recast as a DG method in [9]. It is a Trefftz 
method and its main feature is that all integrals involved in the formulation are boundary integrals, and as such it is much 
cheaper to evaluate than volume integrals from other methods. The idea to couple it with GPWs was proposed in US]. 

4.1 The Ultra-Weak Variational Formulation 

In order to introduce the UWVF, define the slightly more general problem 

f V ■ {AVu) + 2idixdyU - (^1 + ^ + x{x + y)^ u = 0, (U) 

p ■ (AVu) + idiXVyU + uju = Q {—v ■ (TVm) — idiXVyU + lau) + g, (dfl) 


10 













where d = dr ^ idi is a real constant, A = 


dr 

\d\^ 


is therefore a constant matrix, the complex constant fi was defined by 


^ = — 1 — \/l + cr > 0 is a real constant, Q a real valued piecewise constant function on the boundary, and is a Lipschitz 
domain. 

The UWVF is a weak formulation that relies on the use of test functions satisfying the dual equation. Let be a smooth 
test function and be a smooth solution of (31) on an open set O. Then the integration by parts leading to the classical 
weak formulation of (31) yields 


)uv 


/ AVu • Vv — idi j X {dyuv — udyv) + / c(x, y)i 

Jo Jo Jq 

(z/ • {AVu) + tdiXUyu) V = 0. 


(32) 


Now suppose a smooth test function v satisfies the dual equation 


— V • {A\/v) — 2idixdyV + c(x, y)v = 0. 


The equivalent integration by parts yields 


/ AVv • Vu — idi X [dyUv — udyv) + / c(x, y)uv 

Jo Jo Jo 


{u • (AVv) + idiXUyv)u = 0. 


(33) 


(34) 


ao 


Noticing that the volume integral terms in (32) and (34) are identical since A is hermitian, the UWVF contains only boundary 
integrals, thanks to such integration by parts. In order to write explicitly the UWVF for our problem, the following definitions 
are required. 

_ Nh _ 

Consider a Lipschitz domain U C with a mesh U = The boundary of the domain is F = dVt. Let be the 

ki 

diameter of Qk and pk be the maximum of the diameters of the spheres inscribed in O/., and Qk be the value of Q on 
The domain is supposed to be meshed so that Q is constant on each U/.. The mesh is such that 3a such that hk < apk- 
The refinement parameter or mesh size parameter h is then defined hy h = maxh/^. The terminology regular mesh used 
to describe such a mesh comes from [6]. The interface between two mesh elements Qk and Uj, oriented from U/. to Uj, is 
denoted The part of the edge of a mesh element Qk that is part of the boundary of the domain is denoted Tk = FflQU/j;. 

Definition 4.1. On eaeh element Qk of the mesh, the differential operator and its dual are defined on eaeh element of the 
mesh by 

j Lip\Q^ = ^ (-V ■ (AVif)-2tdixdy(p + c(x,y)(p), 

[ ^ (-V ■ (AVip) - 2idixdyr + c{x, y)ip^ , 

and define 

V’) = (fn, AVr-V^- idi fa, x - rdyi>) + fa, c(x, y)r^^ 

(AVtp) + idiXVyr) 

+ idixvyr) 

Ck<p = if 

where a is given by the boundary eondition. 

These definitions of the operators bk and Ck rely on their behavior along interior edges of the domain: 

(Cfel')Sfc,- = 

Denote by (•, •)fe the scalar product on {ip,ip)k = fa, ‘Pi’- Suppose that u is a solution of (31) and u is a solution 
of (33) ; it is now clear that on each element U/. of the mesh one has 

{Lu,v)k — lk{u,v) = — [ —n'{A\/u)v— [ —%diXi'yUV = {bkU,v)k^ 

JdQk 'tdQk 

{u,L*v)k-kiu,v) = - —u-{A'Vv)u+ —idiXUyUV, ( 35 ) 

Jank Jank 


j 

JdQk 


u—u • [AVv) H- idiXUyV = {ckU, bkv)k, 


la 


la 
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which yields, summing over k and using the boundary condition in 

'^{-bku + CkU, blv + Ckv)k - y] '^{-bjU + cju, -blv + Ckv)k 
k k j^k 

-Y.k,Tk^^Qk{-bkU^CkU,-hlv^Ckv)k= ^ {g,-hlv^Ckv)k. 

k,Tk^^ 


(36) 


The function space for the UWVF is V = Y\ while the test function space is defined by 

keli,Nhj 


Nh f 

l-L=Y[^k where Hk = \vk e H^{^k)^ 

k=i ^ 


I/*U/c = 0, 1 

{^k^k T ^k'^k) ^QQj^ ^ L (^dClk) j 


As a result of these definitions, any element of V is actually defined on the edges of every element of the mesh. 


Theorem 4.1. Let u G H^{LL) he a solution of problem ( [M] ) sueh that d^j^u G L‘^{dQk) for any k. Let a >{) he a given 
real number and Q is sueh that Q\dLtk = Qk ^ ^ ond \Qk\ < 1 for all k. Then X G V defined by X\qq^ = Xk with 
Xk (( bk T Ok^u^Qj^f^Q^ satisfies 


T I / -^k{bl +Ck)ek - [ -Xj{-bl+Ck)ek 

k ^ j,j^k ^'^kj ^ ^ 

V / —Xk{-bi^CkYk = E / 


(37) 


fe,rfc5^0‘ 


for any e = (efc)feg|i^jVft] G %■ Conversely, if X is solution of 

of: 

— ^k ^ H (O/0), 
Luk = /|o^, 

{-bk + Ck)uk = Xk, 


then the funetion u defined loeally as the weak solution 


(38) 


is the unique solution of the problem (31). 


This result is classical in the context of UWVF. We refer to O IH |TTJ [121 US]- Equation (37) is the UWVF for (31). 
Note that if u satisfies V • {AVu) G on a Lipschitz domain, then the normal derivative of u on the boundary of that 
domain is well defined. If moreover u satisfies the boundary condition of ( [M| , then the normal derivative of u can be written 
n • Vu = %j\{icFU — idiXUyu) + 0^5 belongs to L^. 


4.2 Discretization 

In order to discretize exactly the UWVF, one would need to use basis functions satisfying the dual equation = 0 but such 


functions are not available. Following the procedure developed in na, we can design GPWs adapted to (pB. These basis 
functions satisfy not exactly the dual equation, but will be designed to satisfy locally the approximation « 0, (f being the 
exponential of a polynomial. More precisely, using Taylor expansions, this subsection focuses on the design of basis functions 
satisfying = 0{h^) for any given order of approximation q. As in [14], the GPWs are defined by the composition of the 
exponential function with a polynomial. As a result, the non-linear system appearing in the design process can be solved for 
the Helmholtz equation. The corresponding algorithm is described here as a more general tool. 


4.2.1 Design of a basis function 

Since the design process is local, focus on a given mesh element Ltk and its centroid G G The question to be answered 
is how to compute the coefficients {Xij) such that the function ip = ^P^^-XG^y-yo) where P = ^ would 

satisfy L'^ip = 0{h^). A natural answer comes from carefully canceling successively the Taylor expansion coefficients of order 
lower than g in !/*(/?. To simplify the notation, define the polynomial 

Vl^ = -dlP - 2drd,,dyP - Idl^d^P - {d^Pf - 2drd^PdyP - \d\HdyPf 

- 2idi{X + XG)dyP 


so that L*ip = Vl* [x - xg, y - yo) + c. 

The system to be solved is then described as follows. The unknowns are the polynomial coefficients (Xij)- for a given deg 
of P, there are P+iKdegP+ 2 ) xhe equations are obtained by equating the Taylor expansion coefficients 
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of !/*(/? of order lower than q to zero: there are equations. Because of the square and product terms, this system 

is not linear and thus requires further investigation. 

Setting degP = q-\-l provides an under-determined system and since then Nu — A^e = 2^^ + 3, 2q^‘i unknowns can be 
fixed conveniently to obtain an invertible system. Since it is not so obvious how to proceed for degP = g, and since the 
system is over-determined for deg P < g, we set from now on deg P = g + 1. The system now reads: for 0 < i j < g — 1, 

0 = —2A2,o ~ 2dr\i^i — 2|(ipAo,2 ~ q ~ 2(i^Ai^oAo,i — 

- 2idiXG\Q,i + c(xG, yc) 

0 = -2A2J - 2dr{j + l)Aij+i - 2idiXG{j + l)Ao,i+i {j > 0) 


^ 1 

- \d\^{j + 2){j + 1)Aoj+ 2 - y] Aij_,Ai,j + -dlc{xG,yG) 

i=o 

j 3 

- 2dr + l)Aij-/Ao,z+i — |dp — k ^ l)Aoj-/c+iAo,/c+i, 

/=0 /c=0 

0 = —(i + 2){i + l)Ai+2,j ~ 2dr{i + l)(g + {i > 0) 

- 2idiXG{j + l)Aij+i — |dp(jf + 2)(jf + l)Aij+2 

- 2idi{j + l)Ai_ij+i + l j^dldlc{xG,yG) 

i 3 

- ^ ^(i - /c + l){k + l)Ai 

k=01=0 

i 3 

- 2dr P P(* - k + 1)(/ + l)Ai -k+l,j-l^k,l+l 

k=0 1=0 
i 3 

1=0 k=0 


(39) 


The next step is to describe which 2g + 3 unknowns can be fixed to make this system conveniently invertible as announced. 
So for a given (i,g), examine the corresponding equation in the system. The key point is to realize that the nonlinear terms 
only involve unknowns Xm,n such that m^n < i + j + 1, whereas the unknowns A^^^ such that m + n = i + g + 2 only appear 
in linear terms. As a result, if for each equation these unknowns involved in the non-linear terms are known, the resulting 
system to be solved would in fact be linear. That is to say the system can be written: for 0 < i + j < g — 1, 


2A2,o + 2drXi^i + 2|(ipAo,2 

= -A^o - 2d^Ai,oAo,i - |dpAg 1 - 2idiXGXo,i +c(xg,^g) 

2A2j- + ‘2‘dr{j + l)Aij+i + |dp(j + 2)(j + l)Aoj-+2 

3 3 

— ‘2iidiXG{^j T l)Ao,ji+i ^ ^ Xi^j—iXi^i 2d'p T l)Aij_/Ao,/-i-i 

1=0 1=0 

^ 1 

- \df ^ + l)Ao,j-/c+lAo,/c+l + vy5^c(xG, ^g), (j > 0) 


/c=0 


{i + 2){i + l)Xi-^2,j + 2dr{i + l)(g + l)Ai+ij+i + \d\‘^{j + 2){j + l)Aij+2 
= -2idiXG{j + l)Aij+i - 2idi{j + l)Ai_ij+i + ^ j^didic{xG,yG) 

i 3 

- ^ ^(i - l){k + l)Ai -k-\-l,j-l^k+l,l 

k=0 1=0 

i 3 

— 2dr P P(i - k + 1)(^ + 1)A, — k-\-l,j—l^k,l-\-l 

k=0 1=0 
i 3 

~ ~ ^ + l)Ai-zj-/c+iA/^/c+i, 

1=0 k=0 


(i > 0). 


(40) 


So first of all, suppose that the right hand side of (40) is known. Now for a given level >C, 0<£<g — 1, there is a 
subsystem of (40) made of the £ + 1 equations corresponding to (i, j) = (i, £ — i). The unknowns of this sub-system are the 
Xi^/:+ 2 -i with 0 < i < £ + 2. So this sub-system has 

• £ + 1 equations. 


• £ + 3 unknowns. 


Moreover the subsystem is tridiagonal since in each equation indexed by i the unknowns involved are A^+ 2 ,£-i, Ai+i,£_i+i 
and Xi^c-i-\- 2 ‘ As a consequence, one easy way to get an invertible subsystem is to fix the two first unknowns Ao ,£+2 and 
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Ai^£+i. The solution of the subsystem then simply reads : for all i from 0 to 


Ai+2,£-i — — 2dr{i + !)(£ — i + l)A4+i 

— \d\‘^{C — i 2){C — i + l)A^^£_i+2) 5 


(41) 


RHS{i) being the right hand side of the corresponding Equation (i,£ — i) in System (40). 

Then suppose the right hand side of a subsystem for a given level C is known. Since this subsystem can be solved thanks 
to ( [4T| , it is clear that the right hand side of the subsystem made of the £ + 2 equations corresponding to (£ j) = (£ £ +1 — i) 
is known. The only remaining step is to initialize this induction process for the level £ = 0. At this level the subsystem is 
only one equation, namely 


2 A 2 ,o + 2 (i^Aiq + 2 |dpAo ,2 = -A^^q “ 2 d^ApoAo,i - \dy>^l^i - 2idiXG\o,i + c{xg, Vg) 


(42) 


So fixing the unknowns Ai^o and Ao,i fixes the right hand side. 

To summarize, for each level £ two unknowns have to be fixed, so 2q unknowns, plus two unknowns for the level £ = 0. 
Notice moreover that the unknown Ao,o does not appear in the system, so it can be fixed without any consequence on the 
previous reasoning. Altogether these are 2^^ + 3 unknowns fixed, corresponding to the difference Nu — between the number 
of unknowns and the number of equations of the initial system. 


Finally, solving system (39) can be described by the following algorithm: 


• fix Ao,o 

• for level £ = 0 


— fix Ai^o and Ao,i 

— compute 

—2dr\i^i — 2|dpAo,2 — A^ 0 “ 2(irAi^oAo,i — 

-2idiXG\{),i +c(xg,^g)) 


V,0 = \ ( 


• for all levels £ from 1 to g — 1 

— fix Ao ,£+2 and 

— compute 

1 


~ 2 y ~ ^ 2)(^ A 1)^0,£+2 “ 2idiXG{C + 1)Ao,£+i 

jC jC 

— 2dr{C + 1)Ai^£+i — Ai^£_/Ai^/ — 2dr ^^(^ + 1)^1,£-Z^0,/+1 

^ 1=0 1=0 

- |d|^ -k + l)Ao,£-fe+iAo,fe+i + ^dyCixG, Vg) ), 

fe =0 ■ / 

Aj+2,£-i = _|_ 2'j(^i -I- 1) ^ ~ “ * + 2)('^ ~ * + l)^i,£-i+2 

(i > 0) —2dr{i + 1)(£ — i ~ 2^d^XG'(£ — i l)\i^c-i+i 

-2idi{C-i + l)Ai_i,£_i+i + t ^ ^^^ dl.dy~^c(xG,yG) 

i C—i 

- y] y](z - k + 1)(A: + l)Ai -/c+l,£-i-zA/c+l,Z 


k=0 1=0 
i C—i 


-2dr yj - k+i){i +i)Ai — k-\-\,C — i —Z'^/c,Z-|-l 

^ k=01=0 

— /c + l){k + l)Xi-i^C-i-k-\-l^l,k-\-l I • 


1=0 k=0 
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4.2.2 Normalization of a set of basis functions 

The choice of polynomial coefficients G {0,1}, 0 < j < g + 1 — i} is then the only remaining step to completely design 

a GPW. This provides a tool to define locally not only one but a set of basis functions. 

First of all, in order to simplify the computational steps, it is natural to set as many coefficients as possible equal to zero. 
Keeping in mind the example of a classical plane wave, the idea proposed here is to fix (Ai^o, Ao,i) - depending on an angle 
0 - to ensure that A 2 ,o = 0? while all the other coefficients of i G {0,1}, 0 < j < g + 1 — i} are set equal to zero. As a 
consequence the resulting GPW ip satisfies p = exp(Ai,o^ + Ao,i^ + O ((x,^)^). From the expression of A 2 ,o provided in the 
algorithm, one can see that this is easily done by fixing Ao,i and setting 


Ai,o = -rfrAo.i + Ag 1 + 2%diXG\o,i) + c{xg, Vg)- 


(43) 


The design of a local set of basis functions is then achieved, following the plane wave example, by fixing the coefficient 
Aq 1 _ ^Jc{xG,yG )_ equi-spaced values of G [0, 27r). 

4.2.3 Summary 

On a given element of the mesh define the basis functions by: 

• k = {xg^Vg) is the center of gravity , 

• p{k) is the number of basis functions, and V/ such that 1 < / < p{k)^ 


- define 9i = - 1), 

— define Ag ^ sindi 


— compute A^ 0 from (43), 

— set the other fixed coefficients 

G {0,1},0 < j < g + 1 - + j 7 ^ 1} to zero, 

• q{k) is the approximation parameter on V/ such that 1 < / < p{k) and \/{i^j) such that 0<i + g <g + l 

— compute the remaining polynomial coefficients according to the algorithm 

— form the corresponding GPW p^^ = exp 

To obtain functions defined on the whole domain ^4, these basis functions {p^]^}i<i<p(^k) are set to be zero on Vt\Vtk. This 
process defines a set of basis functions on f], namely 

£ = Uk£^^{N,p{k),q{k)) where £^^{N,p{k),q{k))P = {<£k }’ 


that will be used to discretize the UWVF 

Thanks to the analysis presented in this section, computing the polynomial coefficients of a GPW does not require to solve 
any system: it reduces to applying the induction formula. All the polynomial coefficients necessary to define the function 
space E can be precomputed and stored, so that the evaluation of any GPW simply requires to read the corresponding 
coefficients from the precomputed table. The implementation that is used to produce the results displayed in the next 
section provides the following timing for p = 7 and g = 4: for a mesh of 34036 elements, the computation of the polynomial 
coefficients table takes 235 seconds while it takes 3 seconds for the classical plane waves equivalent computation ; for a mesh 
of 78574 elements, the computation of the polynomial coefficients table takes 701 seconds while it takes 7 seconds for the 
classical plane waves equivalent computation. These two sizes of mesh are typically the smaller and bigger sizes used in the 
next section. These timings are negligible compared to the time required to build the matrix. 


5 Numerical simulation 

The aim of this section is to show numerical evidence of waves propagating through the mode conversion point. In our model, 
the propagative zones are defined by {x < 0 and y > —x} and {x > 0 and y < —x}. A series of test cases will be designed to 
observe some features of this model, with the concern of observing the incoming wave and avoiding non-physical reflections. 
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Figure 3: Functions /, / and g defined in (44), for i/k = 6. Left: Level curves of /, highlighting the C = /{xk^Vk) and 
level curves together with the line of steepest gradient in the propagative zone. Middle: In the black zone / is constant and 
positive, it is included in the evanescent zone. In the white zone / is constant and negative, it is included in the propagative 
zone. Right: g is non-zero only in the region where \x{x + ^)| < —C. 


5.1 Definition of the test cases 

Since the mode conversion point in the {x^y) coordinate system is the origin, the computational domain is centered at the 
origin. Remember that the medium is propagative \i x{x ^ y) < 0 and evanescent if + ^) > 0. The propagative and 
evanescent zones are separated by two cut-offs, x = 0 and x + ^ = 0. 

As already mentioned, the variable coefficients are key to the mode conversion phenomenon. The two only variable 
coefficients are 


• the zeroth order term, varying as x{x -\-y)^ 

• the first order term, varying as x. 

To observe an incoming wave before and after crossing the mode conversion region, these variable coefficients are set to 
constant values away from the origin. That is to say that the zeroth and first order coefficients will be allowed to vary only 
in a vicinity of the mode conversion point, while they will be set to a constant value in a zone further away from this point. 
In order to introduce only one artificial interface between the constant coefficient zone and the varying coefficient zone, the 
constant zone is chosen to be the same for both coefficients. It is fixed at level curves of the function f{x,y) = x{x -h y), 
for a given transition point {xK^yx)^ see Figure]^ The size of the varying coefficient zone clearly varies with the distance 
between the transition point and the origin. In the propagative zone, that is to say {{x^y) s.t. f{x,y) < 0}, the steepest 
gradient of / is along the line y = —(\/2 l)x. For a given point (xK^yx) along this line and C = f{xK,yx)^ the constant 
zone has four connected components defined by {{x^y) s.t.f{x,y) < C and x < 0}, {{x^y) s.t.f{x,y) < C and x > 0}, 
{{x^y) s.t. f{x,y) > —C and x < 0} and {{x^y) s.t. f{x,y) > —C and x > 0}. In the model studied in this work, the two 
former components are in the propagative zone while the two latter ones are in the evanescent zone. The only remaining 
connected component of the plane is the varying zone. It can also be seen as a transition zone around the mode conversion 
point. In the constant zone the first order term’s coefficient is set to zero, and zeroth order term x{x y) is replaced by 
either ±C in order to be continuous. That is to say that the coefficients used for numerical purposes are 


f{x, y) 


—C if x{x -\-y) > —C, 
Cifx{x^y)< C, g{x,y) 
x{x ^y) otherwise. 


X if \x{x + ^)| < —C, 
0 otherwise. 


(44) 


This function / is continuous while g is discontinuous along the level curves f{x^y) = ±C. Figure shows the modified 
variable coefficients. As a consequence of this setting, a wave propagating from the constant propagative zone toward the 
mode conversion point is expected to propagate until it reaches the varying zone. The wave can then split into a reflected 
wave and a transmitted wave propagating at least partially through the varying zone. Mode conversion then corresponds to 
the existence of a transmitted wave exiting this transition zone on the other connected component of the constant propagative 
zone. 

The antenna is made of a wave guide of width /q and length 4/o, plus a horn, as described in Figure In order to send a 
wave toward the mode conversion point, this antenna is placed facing the origin, in the {x < 0} and {y > —x} region. A given 
transition point (xx^yx) is fixed, defining the amplitude of the transition zone. Then the vertical position of the antenna is 
determined so that if it is along the steepest gradient axis, the distance between the end of the horn and the transition point 
is 12/o- An example of computational domain with the previously described features is also displayed in Figure 

Then the antenna can be placed along any axis with an angle 0 with respect to the y axis, 0 G [0,7r/4], as shown in Figure 
Along the main axis of the antenna, between the antenna and the mode conversion region, the zeroth order coefficient 
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Mode conversion point 




Figure 4: Left: Geometry of the antenna that will be placed in the propagative zone facing the mode conversion point, with 
different incident angles. Right: Computational domain. The blue heavy line represents the boundary of the domain. The 
red dashed lines represent the two cut-off lines, and cross at the mode conversion point. The antenna is placed on the top 
boundary, in the propagative zone. 





Figure 5: Left and middle: Antennas facing the mode conversion point, in the top left quarter of the computational domain. 
The mode conversion point is the bottom right corner of both graphs. The point {xk^ Vk) is the black circle while the mode 
conversion point is red cross. Left: Antenna along the steepest gradient axis. The transition line defined hy x{x ^ y) = 
xk{xk + Vk) and x < 0 is the thick red line. Right: Different positions of the antenna at different angles in the propagative 
zone. Right: Real part of the coefficient c{x{t)^ y{t)) = 1 + l//i + sin^(sin^ — cosO)t^ along the main axis of the antenna, for 
different incident angles, computed with d = 2 — i and = 6. 
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Figure 6: Waves propagating from the antenna placed at an angle 0 = tt/S with the vertical axis, toward the mode conversion 
point (0, 0). The parameter d is d = —2 — i. Left: Computed for i/k = 7. Right: Computed for yx = 12. 


c{x{t), y{t)) = 1 + l//i + /(x(t), ^(t)) is constant on the antenna side and then varying as a quadratic function until the mode 
conversion point. It is odd as a function of t, the parameter describing the position along the steepest gradient line. Figure 
also represents this coefficient c along the main axis of the antenna, for different incident angles 0. The antenna is described 
by x{t) = {sm0)t and y{t) = —(cos6>)t, the parameter t being negative between the antenna and the mode conversion point, 
and 0 at the the mode conversion point. 

The rest of the computational domain is a rectangle aligned with the x and y axis, except for the bottom right corner 
which is cut perpendicularly to the steepest gradient line y = — (v^ + l)x. This cut is meant to avoid artificial reflections of 
waves outgoing from the transition zone toward the {x > 0 and y < —x} region. 

The parameter Q in the boundary condition ( [3l| determines the type of boundary condition. On the bottom wall of 
the antenna, an incoming plane wave, with wave length /q? is sent into the domain by a Robin type boundary condition, 
setting 2 = 0. On the other walls of the antenna, a perfect conductor type boundary condition is set by Q = —1 and g = ^. 
Absorbing boundary conditions are set on on all the other boundaries of the domain, setting Q = 0 and ^ = 0. 

5.2 Results and comments 

A series of results are proposed to illustrate the mode conversion phenomenon in the 2D model introduced in this work. 
Remember that by definition of the O- or X-mode cut-offs, a pure wave of the corresponding type can only propagate on 
one side of the cut-off. As a measure of mode conversion, define the transmission coefficient computed here as the relative 
difference between the maximum wave amplitude before the transition zone, i.e. on {x < 0} D {x{x ^ y) < C}, and after 
the transition zone, i.e. on {x > 0} H {x{x ^ y) < C}. All the following results were computed with p = 7 basis functions 
per element, q = 4 as the approximation order and with either d = —2 — t or d = — tan(7r/8) — i. The absolute value of the 
solution will always be represented with 200 level curves. The transition lines between the constant and varying coefficient 
zones are indicated in blue. Note that the computational domain is meshed automatically by only imposing the mesh size, 
and so the mesh does not resolve the boundaries of the propagative and absorbing region, or the zones of constant and 
variable coefficients. 

To present the type of results obtained. Figure displays two solutions computed with the incident angle 0 = tt/8 and 
d = —2 — 2 , for both yx = 7 and yx = 12. For the yx = 7 case, the antenna parameter is Iq ~ 1.82, the mesh is generated 
for a mesh size of 0.64, and is made of 36548 triangles. For the yx = 12 case, the antenna parameter is /q ~ 1.06, the mesh is 
generated for a mesh size of 0.39, and is made of 63456 triangles. One can observe the same type of behavior: the outgoing 
wave propagates from the antenna, driven by the two cut-offs toward the mode conversion point. Reffection can be observed 
clearly in the transition zone where the coefficients are varying. Comparing these two solutions naturally suggests that the 
size of the transition zone plays an important role in the transmission process. 

The transmission process is expected to depend strongly on the incident angle 0. Figure displays this transmission 
coefficient as a function of the incident angle d, for both d = —2 — i and d = —tan(7r/8) — %. Typical values of various 
geometry and mesh parameters used to produce Figure are given in Table As expected, the transmission is very low 
when d is close to zero and 7r/4, which correspond to launching the wave parallel one of the cut-offs, while a peak is observed 
between those two regimes, which corresponds to the wave propagation direction closer to the steepest gradient line of the 
variable coefficients. Figure [^displays the solution with the maximum transmission as well as a solution completely reflected 
by the cut-off before reaching the mode conversion region. 

Figure [^evidences the importance of the size of the transition zone on the mode conversion efficiency. Indeed the maximum 
transmission coefficient is a decreasing function of yxi which is a measure of the size of the transition zone. Moreover the 
size of the transition zone also influences the angle at which this maximum coefficient occurs. See Tables and 
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Figure 7: Transmission coefficient as a function of the incident angle 0^ for different sizes of transition zones. Left: Computed 
for d = —2 — i. Right: Computed for d = — tan(37r/8) — i. 


Table 1: Geometry and mesh parameters for different sizes of the transition zone. 


Vk 

6 

7 

8 

9 

10 

11 

lo 

2.13 

1.82 

1.59 

1.42 

1.28 

1.16 

mesh size 

0.73 

0.64 

0.57 

0.46 

0.43 

0.39 

number of triangles 

25157 

36548 

42412 

60944 

66511 

78470 


Acknowledgment 

Many thanks to Harold Weitzner for introducing me to mode conversion and guiding me along the way, and to Jonathan 
Goodman for our many discussions. I also thank Teemu Luostari for providing his 2D PW-UWVF code for elasticity 
equations. 

This work was supported in part by the U.S. Department of Energy, Office of Sci- ence, under Grant No. DE-FG02- 
86ER53223. 

References 

[1] I. Babuska and S.A. Sauter, Is the pollution effect of the EEM avoidable for the Helmholtz’s equation considering high 
wave numbers? SIAM Rev. 42 (2000) 451-484 

[2] P.T. Bonoli Review of recent experimental and modeling progress in the lower hybrid range of frequencies at ITER 
relevant parameters AIP Conference Proceedings 1580 15-24 http://dx.doi.Org/10.1063/l.4864497 

[3] P.T. Bonoli Electromagnetic mode conversion: understanding waves that suddenly change their nature Journal of 
Physics: Conference Series 16 (2005) 35-39 

[4] A. Buffa and P. Monk, Error estimates for the Ultra Weak Variational Eormulation of the Helmholtz equation, ESAIM: 
Mathematical Modelling and Numerical Analysis42 (2008) 925-940 




Eigure 8: Comparison of two solutions, highlighting the influence of the incident angle on the transmission coefficient. Left: 
Computed for yx = 10, d = —2 — 2 and 0 = 0.74, corresponding to a transmission coefficient T ^ 10“^. Right: Computed 
for = 6, d = — tan(37r/8) — i and 0 = 0.371, corresponding to the transmission coefficient T = 0.065. 


19 
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